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ABSTRAGT 

Extremely large surveys with future experiments like Euclid and the SKA will soon allow us to 
access perturbation modes close to the Hubble scale, with wavenumbers k ^ H. If a modified gravity 
theory is responsible for cosmic acceleration, the Hubble scale is a natural regime for deviations from 
General Relativity (GR) to become manifest. The majority of studies to date have concentrated on the 
consequences of alternative gravity theories for the subhorizon, quasi-static regime, however. In this 
paper we investigate how modifications to the gravitational field equations affect perturbations around 
the Hubble scale. We choose functional forms to represent the generic scale-dependent behaviour of 
gravity theories that modify GR at long wavelengths, and study the resulting deviations of ultra 
large-scale relativistic observables from their GR behaviour. We find that these are small unless 
modifications to the field equations are drastic. The angular dependence and redshift evolution of the 
deviations is highly parameterisation- and survey-dependent, however, and so they are possibly a rich 
source of modified gravity phenomenology if they can be measured. 


1. INTRODUGTION 

The past decade of work on extensions of General Rel¬ 
ativity (GR) has primarily been motivated by the hope 
of explaining cosmic acceleration without a fine-tuned 
cosmological constant, A. The inferred energy scale of 
A is approximately the Hubble scale today, fc ~ iJo ~ 
2 X 10““^ Mpc“^; so, for theories hoping to dispose of 
A, it seems natural that the cosmological horizon should 
act as a threshold at which corrections to GR become 
relevant. Studying how modified gravity (MG) theories 
affect Hubble-scale perturbation modes could therefore 
offer interesting insights into the dark energy and cos¬ 
mological constant problems. 

Long wavelength (infrared; IR) modifications of grav¬ 
ity have mostly been explored in the context of 
braneworld models (Dvali et al. 2000; Maartens 2004), 
non-local theories (Maggiore & Mancarella 2014; Dirian 
et al. 2014; Gonroy et al. 2015) and massive grav¬ 
ity/bigravity (de Rham et al. 2011; D’Amico et al. 2011; 
Hassan & Rosen 2012). In some theories the features 
causing the long-wavelength modifications also result 
in ultraviolet (UV) phenomenology, such as the mod¬ 
ification of forces at a source-dependent distance r* 
(Babichev & Deffayet 2013), extra graviton polarisations 
(Dvali 2004), and the existence of a UV cut-off to the 
theory. In this paper we will focus solely on the IR 
aspects of modifications to gravity; any associated UV 
phenomenology is unlikely to affect the observables we 
consider here. 

In contrast to the demanding simulation requirements 
and astrophysical uncertainties that plague non-linear 
wave modes k > 0.1 Mpc“^, ultra-large scales are a rel¬ 
atively ‘clean’ and safely linear regime in which to test 
gravity. Perturbation modes with wavelengths approach- 
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ing the Hubble scale are mostly beyond the reach of cur¬ 
rent large-scale structure experiments however (although 
see, e.g. Leistedt et al. 2013; Hernandez-Monteagudo 
et al. 2014), and so the bulk of work on constraining al¬ 
ternative gravity theories has focussed on significantly 
subhorizon scales, often using quasi-static approxima¬ 
tions to simplify calculations (Amendola et al. 2013b; 
Silvestri et al. 2013; Sawicki & Bellini 2015). The obser¬ 
vational situation is soon set to change though, as forth¬ 
coming surveys of large-scale structure - notably wide- 
angle galaxy and 21cm intensity mapping surveys (LSST 
Science Collaboration 2009; Schlegel et al. 2011; Amen¬ 
dola et al. 2013a; Bull et al. 2015; Jarvis et al. 2015) - 
will make it possible to probe ultra-large scales across gi¬ 
gantic survey volumes at z < 5 (Yoo et al. 2012; Camera 
et al. 2013; Camera et al. 2015; Raccanelli et al. 2015; 
Alonso et al. 2015). 

The purpose of this article is to understand how 
Hubble-scale gravitational corrections affect various ob¬ 
servables on ultra-large scales. The space of MG theo¬ 
ries is complex, and there is no consensus on a favoured 
model. Hence, as far as possible, we will draw only 
upon generic physical features (Lorentz invariance, mass 
scales, etc.) instead of adopting a particular theory. 
From §3.4 onwards we choose some specific ansatzes in 
order to obtain example numerical results; our calcula¬ 
tions up to this point are general, however. 

Horizon-scale modes require careful handling: cal¬ 
culations performed in different gauges will no longer 
agree in general^, and relativistic effects become impor¬ 
tant. The relativistic corrections impacting horizon-scale 
large-scale structure observables in GR were worked out 
in Yoo et al. (2009); Yoo (2010); Bonvin & Durrer (2011); 
Challinor & Lewis (2011); Yoo et al. (2012); Jeong et al. 
(2012); in §4 we will consider how they are affected by 

^ Unless the quantity being calculated is a true observable. 
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Fig. 1. — Comoving Hubble scale, = "H, as a function of 
redshift and comoving distance to the last scattering surface (blue 
solid/dashed lines). Also shown are the (Limber-approximated) 
comoving scales probed by a given multipole, £ ~ /cr (black 
solid/dotted lines), the matter-radiation equality scale, k^q (yellow 
line), and the normalised fractional energy density in dark energy, 
^DE(-2)/f^DE,o (horizontal red lines), (c.f. Zibin & Moss 2014) 


modifications to the gravitational field equations. De- 
tection of the relativistic effects will be challenging (Yoo 
et al. 2012; Maartens et al. 2013; Alonso et al. 2015), as 
their impact is limited to low multipoles of the angular 
power spectrum of the galaxy distribution (see Fig. 1). 
A number of other physical and systematic effects can 
also contaminate the signal (Bruni et al. 2012; Pullen & 
Hirata 2013; Leistedt et al. 2013; Alonso et al. 2015). 
Nevertheless, it is important to take the relativistic cor¬ 
rections into account when studying the Hubble scale. 

This paper is structured as follows. After motivating 
our effective description of the gravitational field equa¬ 
tions (§ 2 ), we first build intuition about how the dark 
matter density perturbation is altered by horizon-scale 
deviations from GR (§3). We then apply these ideas to a 
large-scale structure observable, the angular power spec¬ 
trum of source number counts (§4), adopting a formula¬ 
tion that can describe both number counts and bright¬ 
ness temperature fluctuations simultaneously. We dis¬ 
cuss our results and their observational feasibility in §5. 

2 . PARAMETERISATION 

Observations to date are consistent with a constant 
dark energy equation of state of ru = —1 (Planck Col¬ 
laboration XIV 2015), leaving the evolution of perturba¬ 
tions as the most likely arena for seeing deviations from 
ACDM on cosmological scales. As such, we will assume 
a cosmological expansion history that is indistinguish¬ 
able from ACDM in what follows. This restriction of a 
fixed background can easily be relaxed at the expense of 
introducing more parameters into our framework. 

The ultra-large scale modes that we are interested in 
are very well-described by linear perturbation theory. 
Studies focussing on the quasi-static regime have shown 
that gravity theories endowed with a single scalar degree 
of freedom are encapsulated by the following effective de¬ 
scription on small (but still linear) scales (see §3.1 for our 


conventions regarding the metric potentials):^ 

— 2 fc^T = 87 rGAra^pD/i (a, fc) ( 1 ) 

$ = 4'7(a,fc), ( 2 ) 

where D is a gauge-invariant fractional density perturba¬ 
tion, and a sum over fluid species is implied on the RHS 
of Eq. (1). (Note that D here is often denoted by A in 
other works; we follow the conventions of Durrer (2005) 
used by the CLASS code in §4). The hats on the metric 
potentials indicate fully gauge-invariant quantities (the 
Bardeen potentials), which we introduce for later conve¬ 
nience. {( 1 , 7 } are two functions of conformal time and 
scale. For theories with second-order equations of mo¬ 
tion (e.o.m.s), their ^-dependence is known; it has been 
derived for particular models in de Felice et al. (2011); 
Gleyzes et al. (2013); Bloomfield (2013), but is in fact a 
straightforward consequence of a few general principles: 
gauge invariance, Lorentz invariance, and second order 
e.o.m.s. (Baker et al. 2014). 

With a little care, the same arguments can be used 
to justify the use of Eqs. (1) and (2) on scales beyond 
the quasi-static regime too. The scale-dependence of jl 
and 7 is a little more complex here, as terms that go like 
Ti/k cannot be neglected. Similarly, one must allow for 
the existence of a new mass scale (M) and/or timescale 
(F) that characterise deviations from GR, and these will 
enter the form of {/i, 7 }.^ 

The arguments presented in Baker et al. (2014) moti¬ 
vate the following approximate forms for use outside the 
quasi-static regime: 

_ 0(1) + o(U)+o(e;)+o(^)+... 

o(i) + o(a;) + o(g) + o(¥) + ... 

( 3 ) 

It should be understood that all terms are multiplied by 
time-dependent coefficients specific to a particular grav¬ 
ity theory; in model-independent calculations, a sensible 
ansatz must be supplied for these. The key point is that 
fl and 7 have the same fc-dependence as their quasi-static 
counterparts, with some additional scales F and M com¬ 
ing into play. Given that the purpose of this paper is 
to investigate the effects of Hubble-scale corrections to 
GR, we will take F ~ M ~ H, which further simplifies 
Eq. (3). 

This line of reasoning justifies the use of a simple 
parameterisation such as Eqs. (1) and (2) to explore 
large-scale MG phenomenology, which we treat in the 
next section. An earlier model-independent treatment 
of horizon scales was incorporated in the Parametrized 
Post-Friedmann formalism of Hu & Sawicki (2007), but 
the parameterisation we use here has the advantages 
of simplicity and updated knowledge about the generic 
scale-dependence of MG theories, as discussed above. A 
scheme for extending beyond the quasi-static regime in 
Horndeski gravity was very recently presented in Lom- 

^ Eq. (1) is not strictly the Poisson equation, as it has the time¬ 
like gravitational potential on the LHS (i.e. if instead of $). 

^ The usual quasi-static approximation assumes that T k, 
allowing such time derivative terms to be neglected; on larger, non¬ 
quasistatic scales we must keep them. 
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briser & Taylor (2015), using ideas similar to those dis¬ 
cussed above. ^ 

In §3.4 we will evaluate our results numerically, using 
a simple model based on Eq. (3). We do not claim that 
Eqs. (1-3) apply comprehensively to all gravity theories 
on large scales,^ only that they are sufficiently general for 
use in a model-independent investigation such as ours. 
Similar forms were used (also outside the quasi-static 
regime) in recent constraints on dark energy by Planck 
Collaboration XIV (2015), for example. Constraints on 
{/i, 7 } in two bins of k, using data from the WiggleZ, 
BOSS, and 6 dE surveys, were recently presented in John¬ 
son et al. (2015). 


3. SCALE-DEPENDENT GROWTH 

Even with relativistic corrections taken into account, 
the matter density perturbation Dm remains the dom¬ 
inant contribution to the observables considered in §4. 
In this section, we will build some basic analytic intu¬ 
ition about the effects of our modified equations ( 1 ) and 
(2) on the evolution and scale-dependence of Dm- (We 
note that the growth rate of matter density perturbations 
on large scales in a few specific gravity models and the 
Parametrized Post-Priedmann framework was considered 
in Lombriser et al. 2013). 

Issues of gauge choice become crucial when dealing 
with non-observable quantities on near-horizon scales, so 
we will work with gauge-invariant variables as far as pos¬ 
sible. To avoid cluttered expressions, we will suppress the 
arguments (fc, a) for most variables. 


3.1. Set-up 

Let the line element showing the perturbations of the 
FLRW metric in a general gauge be (Durrer 2005; Bonvin 
& Durrer 2011): 

= 0 ^( 77 ) I — (1 -I- 2A)dr]^ — 2Bidx'^dr] 

+ [(1 + 2iLi)Jy + 2HTrj] dx^dx^y (4) 

As we are interested in only the spin-0 components 
of the perturbations here, we can write = —kiB/\k\ 

and HtiJ = {kikj — ^5ij)HT/k"^■ For a fluid with four- 
velocity we define the spin -0 component of its spatial 
3-velocity perturbation to be: 


u = 


a 



(5) 


With these definitions in hand, we will work in terms 


The authors use the evolution at a fixed ‘pivot scale’, fc*, to 
define a timescale equivalent to our F; see Eqs. (3.1) and (3.2) of 
Lombriser & Taylor (2015). 

® Their derivation relies on assumptions such as that the evolu¬ 
tion of a scalar degree of freedom can approximately be described 
by X ~ TxX- 


of the following gauge-invariant variables: 


4 - = 

$ = 

D = 

V = 

Y = 





S -h 3(1 -h w)^(v - B) 
k 


Hrj 


V — 



D-3il + w)\^ 
D -\- 3(1 -I- w)TZ 




( 6 ) 

(7) 

( 8 ) 

(9) 

( 10 ) 

( 11 ) 


where the last line above defines Tico, which becomes 
equal to (minus) the curvature perturbation in the co¬ 
moving gauge. All of our calculations will concern times 
well after matter-radiation equality, such that pressure¬ 
less matter (indicated by a subscript M) is the domi¬ 
nant clustering matter component and anisotropic stress 
is negligible. 

We need four equations to solve for the four indepen¬ 
dent variables $, 4^, Dm and Vm (from Eq. (10) we see 
that 1 m is simply a linear combination of these). Two are 
provided by Eqs. (1) and (2), where it should be noted 
that {/i, 7 } are defined in a gauge-independent manner. 
The other two are the standard fluid equations for pres¬ 
sureless matter perturbations: 


Ym = -kVM I4i = -nVM + ki>, (12) 

where overdots denote conformal time derivatives. The 
usual Jm term is sitting inside our variable Ym- Note that 
we have assumed here that pressureless matter follows its 
standard conservation laws, thereby excluding some cou¬ 
pled dark energy models (for example) from our treat¬ 
ment. We leave such extensions to a future work; note 
that there will be some resulting changes to Eq. (28). 


3.2. Evolution o/Ym and Dm 

The ‘cleanest’ variable to solve for is Ym. If we solve 
for Ym instead of Dm directly, we avoid the need to take 
time derivatives of our parameterisation functions fi and 
7 ; this is advantageous both for simplicity’s sake, and 
also because p, and 7 are effective, rather than exact, 
objects. 

Substituting Eqs. (2), (10) and the second of Eqs. (12) 
into Eq. (1), and using 3 'H^Dm = ^t^Gno^pm, we obtain 

(Ym - 3^Ym) 

-2kH = 3n^nMP > 3^0 (13) 

(1 + ^^VImPi) 

Note that the matter fraction here is time-dependent; we 
will use Dmo for its value today. Comparing Eqs. (13) 
and ( 1 ), we infer that 
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Differentiating the first of Eqs. (12) and using Eqs. (2), 
(13) and the second of (12) then feads us to an evolution 
equation for Ym, 


+ HYm 


1 + 


3 

2 


fi-Yu 


1 

1 + 


= 0 . 


(15) 


The solutions of this equation (as a function of k and z 
or a) can then be used to find the evolution of the matter 
density perturbation, Dm, using Eq. (14). 


3.3. Sub- and Super-Horizon Limits 

Using the parameterised Poisson equation, Eq. (10) 
can be rewritten as 


1 m = -Dm 1 + 




2 fc 2 J 


-SfVu. (16) 


In the subhorizon limit, k H, this gives Fm —t -Dm,® 
which is a manifestation of the fact that gauge differ¬ 
ences are generally unimportant on small distance scales. 
In this limit, Eq. (15) recovers the usual evolution equa¬ 
tion for Dm on small scales (e.g. Amendola & Tsujikawa 
2010; Baker et al. 2014), 

iii + HlM-^H^D m/IFm = 0. (17) 

We see that small scales are sensitive only to fl and not 
7 , with enhanced growth (relative to GR) expected for 
values /i > 1 . 

To understand the behaviour of superhorizon scales, 
let us first consider the GR case. In the limit k and 

= 7 = 1 , Eq. (15) becomes 

Y°^ + 2nY°^ ~ 0 , (18) 

that is, the equation for Fm becomes unsourced and 
hence is solved by Fj^^ —>• constant on ultra-large scales; 
see Fig. 2. (The second solution for the above equation 
blows up as a —)■ 0 , and so must be removed by choice of 
boundary conditions.) 

To proceed, we now consider small deviations of quan¬ 
tities from their GR values, writing Fm = Fm^ -|- JFm, 
Dm = -Dm^ -I- i5Dm, ^ = 1 -\- Sp and 7 = 1 -I- Jy. For 
brevity we will not display the full linearised versions of 
Eqs. (14) and (15), but their large-scale limits are 


5Fm + 2nSYM - HYonSl oz 0 
6 Dm (5Fm 


r,GR 


vGR 

-^M 


- (J/i-bJy). 


(19) 

( 20 ) 


Note that although Fj^^ —>■ 0 for fc —)■ 0, the first term 
in Eq. (20) does not have to diverge if <51 m is decreasing 
with fc at a similar rate. In Appendix A, we show that 
this ratio is given on ultra-large scales by 


J {n6j)d7^. ( 21 ) 

-"M 


® Assuming that /iy becomes a scale-independent function of 
time in this limit, which is true of Eqs. (3), (22) and (23). 



10'' 10'^ 10'^ 10'^ 10° 

a 


Fig. 2.— Growth of the variable Ym (defined in Eq. 10) in GR, for 
different wavemodes. Note that the large-scale mode k = 2 X 10“^ 
Mpc“^ scarcely evolves. Curves have been normalised to their 
initial values, which depends on k. 


From Eq. (20) we see that the large-scale behaviour of 
the dark matter density perturbation Dm is controlled by 
the interplay of two terms. For a late-time modification 
of gravity, one expects the integral in Eq. (21) to be the 
subdominant term, as the integrand is suppressed by the 
decreasing Hubble factor just when Jy begins to grow. 
These formulae will help us to understand the large-scale 
behaviour of the two numerical examples presented in the 
next subsection. 

Note that whilst small scales were only sensitive to 
/i, large scales depend on both /I and 7 via Eq. (20). 
Also note that the influence of p changes directionality: 
on small scales, p, > 1 enhances the growth of GDM 
perturbations, whilst on ultra-large scales it suppresses 
growth relative to the GR case (this feature was also 
noticed in Hall et al. 2013). 

In between the the large- and small- scale limits exists 
a transition regime, requiring the full numerical solution 
of Eq. (15). It is this regime which is the most interesting 
observationally - even the next generation of surveys will 
not probe scales large enough to see the k limit at 
late times (corresponding to the area left of the solid blue 
line in Fig. 1). In the next subsection we solve for the 
full range of fc for two contrasting examples of {/ 1 , 7 }. 

3.4. Example Ansatz 

The derivations of the previous section are independent 
of the precise forms of the parameterisation functions p 
and 7 . To solve for the full near-horizon behaviour, we 
now need to pick concrete forms for p and 7 . 

As argued in §2 and Baker et al. (2014), a motivated 
scale-dependent form of p and 7 (for theories with one 
new d.o.f.) is a ratio of polynomials in {M/k^ and 
(r/fc)^, where M and T are a mass scale and timescale 
characteristic of any modifications to GR. In this paper 
we wish to explore scenarios that preserve GR at short 
distances whilst allowing new phenomena to emerge in 
the regime k ^ LL. For the present case, then, we simply 
take M T ^ H. 

The precise form of Eq. (3) is difficult to explore in a 
model-independent manner, as one then has no guidance 
for how the time-dependent coefficients are related to one 
another. We will therefore use a slightly altered version 
of this scale-dependence that makes preservation of the 
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Fig. 3. — Growth of Ym and Dm? as defined in Eqs. ( 8 ) and (10). The left-hand panel employs the late-time emergence (A-type) ansatz, 
the right-hand panel uses the early-time ansatz (‘CPL-like’), see Eqs. (24-26). The dark and pale curves in the two lower panels correspond 
to cases with Sfi and S'y switched on individually, respectively. The top panel shows the time- and scale- dependence of either fi or 7 , which 
control the deviations from GR. Note that although in the CPL-like case jl{z = 0) = 1, 7(2 = 0) = 1, there is an integrated effect from the 
modified field equations that persists to late times. 


GR limit more transparent (a similar form was used in 
Eqs. (46) and (47) of Planck Collaboration XIV 2015): 


/i(a, fc) 


7(a,fc) 


1 + //i(«) 


1 + fj{a) 




( 22 ) 

(23) 


Here, fp,{a) and 7(a) are functions of scale factor that 
control the time evolution of MG effects, and A is a pa¬ 
rameter that determines the location (in Fourier space) 
of the transition from the small-scale GR regime to the 
large-scale MG regime. We will adopt A = 10 in this 
work, as this allows MG effects to become important 
only on scales larger than fceq, which is effectively the 
largest scale probed by large-scale structure surveys to 
date. While picking A = 1 may seem more natural, this 
would limit deviations from GR to multipoles £ < 3, as 
shown in Fig. 1. (We note in passing that subhorizon 
transition scales can occur in models of a clustering dark 
energy component with subluminal sound speed.) 


We must also choose some parameterisations of the 
time-dependence of 7(a) and 7 (a) to explore. We com¬ 
pare two ansatzes that have opposite behaviours; one 
where the MG functions grow larger with time, and an¬ 
other in which they decay. This will allow us to ex¬ 
plore the interplay between different terms on ultra-large 
scales found in Eq. (20). To facilitate comparison with 
the existing literature, we used the two time-dependent 
ansatzes considered by Planck Gollaboration XIV (2015) 
for their constraints on dark energy and gravity: 

1. A-like/late-time emergence (recall that we invoke a 
ACDM expansion history): 

/A(a) = /AA%k), = (24) 

“AO “AO 

2. Behaviour analogous to the Chevallier-Polarski- 
Linder (GPL; Chevallier & Polarski 2001; Linder 
2003) equation of state parameterisation: 

//i(«) =//io +//ii(l - a) (25) 

/7(«) =/70 +/7i(l - a)- (26) 
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Option (1) is the natural choice for MG effects that hope 
to explain late-time acceleration, whilst option ( 2 ) hy¬ 
pothesises that any deviation from an equation of state 
of w = — 1 might be accompanied by modihed perturba¬ 
tions. Given the tight constraints on early dark energy 
from the GMB (Galabrese et al. 2011), we impose an 
additional cutoff in our numerical implementation that 
smoothly brings the modifications to zero at z > 5. 

Fig. 3 shows the deviations of Ym and Dm as ratios of 
their GR values for both ansatzes, as a function of k at 
z = 0. For the A-like ansatz (left panel), two cases are 
plotted: {//iA, /ja} = {0-1,0}, {0,0.1}. We have chosen 
to vary fi and 7 individually in order show their distinct 
effects; in a realistic MG model one would typically ex¬ 
pect them to vary simultaneously. The 10% amplitude 
of deviations from GR was chosen simply for illustration. 

The right panel side of Fig. 3 shows analogous plots for 
the GPL-like ansatz. The curves plotted correspond to 
{hi, hi} = { 0 . 1 , 0 }, { 0 , 0 . 1 }, with {/^o, /^o} = ,{ 0 , 0 } 
for both; once again, this particular example is motivated 
only by simplicity. 

The top row in both panels shows the time- and scale- 
dependent form adopted by either fl or 7 . Note that the 
transition away from the GR regime extends to wave 
numbers some way below the horizon. The following 
other features are also of interest: 

• The ultra large-scale behaviour of the matter den¬ 
sity perturbation can be understood from Eqs. (20) 
and (21). In Fig. 3, where the MG effects be¬ 
come relevant at late times, the integral term in 
Eq. (20) is subdominant, so the large-scale value of 
Dm/D^ is primarily set by the maximum of 5jl 
or 5^. (The offset between the curves is due to the 
integral contribution to 7 .) 

• In the GPL-like case, the integrand of Eq. (21) is 
larger, and this dominates the large-scale deviation 
if 5^ 7 ^ 0. Since i5/i —>■ 0 as z —)■ 0, the deviations 
vanish on large scales for {//ii, fji} = { 0 . 1 , 0 }. 

• As mentioned in §3.3, horizon-scale modifications 
of the effective gravitational constant suppress 
growth; this contrasts with the enhanced growth 
they cause on subhorizon scales. 

• The time-dependence of MG effects plays a strong 
role in their impact on the density perturbations; 
compare the ~ 10% deviations for the A-like ansatz 
with only ~ 2% deviations for the GPL one. 

Much of this (qualitative) behaviour is mirrored in the 
observable angular power spectra calculated in §4, as the 
density-density (and density-RSD) term is dominant over 
a wide range of scales and redshifts (Bonvin & Durrer 
2011 ). 


3.5. Comoving Curvature Perturbation 
Starting from Eqs. (1), (2) and (12), one can show that 
the comoving curvature perturbation TZco = — VH/k 

is conserved on ultra-large scales, provided that jl does 
not tend to zero there. We relegate the derivation to 
Appendix B. This is consistent with the well-known re¬ 
sults of Wands et al. (2000) and Bertschinger (2006) , who 


showed that this feature should hold true in any metric 
theory of gravity. 

The conservation of TZco, together with the conserva¬ 
tion equations for CDM, immediately imply that the evo¬ 
lution of the metric potentials on superhorizon scales is 
given by 




( mu 

\ n‘^-h ) 


^^ {nh-2ih) 

{W-H) 


= 0 . 
(27) 


As noted in Bertschinger (2006) and Hall et al. (2013), 
this equation completely fixes the evolution of superhori¬ 
zon potentials in terms of the cosmological expansion 
rate once a relationship between <!' and 'L has been spec¬ 
ified. If 7 7 ^ 1, this means that deviations from GR are 
expected on superhorizon scales, even if the background 
expansion history is identical to that of AGDM. 

Notice that it is the slip ratio, 7 , that is most relevant 
to superhorizon behaviour (see also Eq. 19), in contrast 
to the controlling influence of jl on small scales (Eq. 17). 


4. ULTRA-LARGE SGALE OBSERVABLES 

Observables such as galaxy number counts are primar¬ 
ily sensitive to fluctuations in matter density, which we 
explored in §3, but a number of other effects can also 
contribute corrections to the observed quantity. Lensing 
by intervening matter, peculiar velocities of the sources, 
evolution of potentials along the line of sight, and sev¬ 
eral purely relativistic effects lead to distortions of the 
redshift-space positions of objects, their flux distribu¬ 
tion, and the survey volume. All of these effects must be 
taken into account in a fully correct treatment - espe¬ 
cially on ultra-large scales, where the corrections can be 
large. Complete calculations up to first and even second 
order in perturbations have been performed for several 
observables (Bonvin et al. 2006; Bonvin & Durrer 2011; 
Challinor & Lewis 2011; Jeong et al. 2012; Yoo et al. 
2012; Umeh et al. 2012; Hall et al. 2013; Andrianomena 
et al. 2014), some of which have been implemented in 
the CLASS and CAMB codes (Challinor & Lewis 2011; 
Di Dio et al. 2013). 

Previous work on tests of GR on ultra large-scales has 
mostly focused on the integrated Sachs-Wolfe (ISW) ef¬ 
fect (e.g. Lue et al. 2004; Zhang 2006; Song et al. 2007), 
which is a sub-dominant contribution to the relativistic 
observables. More extensive treatments have been given 
by Yoo et al. (2012); Lombriser et al. (2013) and Hall 
et al. (2013), who calculate the full linear relativistic ex¬ 
pressions for perturbed source number counts and 21 cm 
brightness temperature fluctuations respectively. The 
former propose to test GR by separately measuring the 
relativistic ‘correction’ terms that are typically neglected 
on subhorizon scales, and calculate these effects for sev¬ 
eral different MG theories in Lombriser et al. (2013). 
Hall et al. (2013) specialise to f{R) gravity, and per¬ 
form a principal component analysis of the 21 cm power 
spectrum (although they do not specihcally separate out 
large-scale modes). 

We will take a different approach, focusing exclusively 
on Hubble-scale corrections to GR using the model- 
independent ansatzes of Eqs. (24-26). The observable 
that we consider is the total linear perturbation to the 
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Fig. 4.— Fractional deviation of the number count power spectrum from its GR behaviour, for the A (left) and CPL-like (right) 
parameterisations. Four redshift bins are plotted for each set of parameters, all with width A 2 : = 0.1. We have fixed A = 10. 


number counts of sources in a flux-limited redshift survey 
(Di Dio et al. 2013), 


AN(ft,fc,z) = &DM + -nfVM (28) 


+ ^' + (5s-2)$ + 'H"^4> 
2-5s 
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rrs 

/ dr 

rq — r 
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($ + 5*) 
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■H2 rsn 


+ 5s - 


X ( 'I' + i(k • n)VM + / dr ($ + 4') I . 


Here, b is the linear bias, s is the magnification bias (the 
logarithmic derivative of the source number density with 
respect to limiting magnitude), is the evolution bias 
(which describes the proper evolution of the source num¬ 
ber density), and Vn is the angular Laplacian. Subscript 
S denotes evaluation at the source. The expression for 
the fractional perturbation to the 21 cm brightness tem¬ 
perature can be obtained directly from this expression 
simply by setting s = 2/5 (Hall et al. 2013). Also see 
Bruni et al. (2012) for a discussion of how to define the 
bias on ultra-large scales. 

Many of the terms in Eq. (28) are well known, and 
some have already been measured by large-scale struc¬ 
ture surveys. The first line consists of the source density 
(including a gauge correction) and a redshift-space dis¬ 
tortion term (oc dr(V ■ n)). The whole third line is the 
lensing convergence, and the final term on the last line 
is the integrated Sachs-Wolfe (ISW) effect, oc ($ - 1 - 4 '). 

Most of the other terms are strongly suppressed in the 
subhorizon, quasi-static regime, as they scale as H/k or 


("H/fc)^ with respect to the dominant density term. As 
such, they can be neglected to a very good approximation 
by essentially all large-scale structure surveys to date. 
As discussed previously, this will no longer be the case 
with future surveys however, which will be able to access 
considerably larger scales. 

Each of the relativistic correction terms will be modi¬ 
fied in different ways by the introduction of modifications 
to GR, and so the sum of these terms could potentially 
exhibit interesting features in a MG theory. The ques¬ 
tion of how the total observable is modified, and whether 
those modifications can be measured, is the focus of the 
rest of this section. 

4.1. Calculating relativistic observables in MG 

We now wish to obtain the full linear relativistic ex¬ 
pression for the number counts, Eq. (28), in the presence 
of modifications to GR represented by our Eqs. (1-3). 
The derivations in Bonvin & Durrer (2011) and Ghalli- 
nor & Lewis (2011) did not assume the GR field equa¬ 
tions, and used only perturbations to purely geometric 
quantities like the Jacobian and luminosity distance, plus 
stress-energy conservation. As such, the MG parameter- 
isation described above affects An only by changing the 
source terms in Eq. (28), e.g. by modifying the solutions 
for ^{k,rj) and 4>(fc, 77 ), leaving the form of the equa¬ 
tion itself unchanged. This is not completely general - 
for example, the assumption of normal stress-energy con¬ 
servation for matter excludes some gravity theories, like 
those with non-universal matter couplings (Amendola 
2000; Pettorino & Baccigalnpi 2008; Brax et al. 2010). 
It is precisely these theories for which the parameterisa- 
tion of Eqs. (1) and (2) is ill-suited though, and so we 
will not consider them here. Note that a more general 
expression, valid for any metric theory of gravity, can be 
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Fig. 5.— Left panel: Angular power spectra for individual terms from Eq. (28), for GR, relative to the total number count power spectrum 
(see right panels for key). Right panels: Fractional deviation of the number count power spectrum for the A parameterisation, as in Fig. 4, 
for the same terms from Eq. (28). 


obtained with only relatively minor changes to Eq. (28) 
(Bonvin 2014). 

Up to the assumptions we have made, then, all that 
is needed to calculate the number count perturbation 
in MG theories is to modify the GR calculation of 
Eq. (28) to use the MG solutions for the source terms 
$, th, Umj and We implement this by modifying 

the GLASS Boltzmann code^ to calculate the angular 
power spectrum of An using the MG perturbation equa¬ 
tions, Eqs. (1), (2), and (Dl). We will consider only late 
time modifications to GR, which allows us to take GR 
initial conditions at high redshift, and to leave the GMB 
calculation unmodified. Details of our implementation 
are given in Appendix D. 

4.2. MG effects on ultra-large scale observables 

Eig. 4 shows the fractional difference of the number 
count angular power spectrum, Ce{z,z), from GR, for 
four different choices of MG ansatz/parameters. The 
galaxy bias was taken to be constant, 6 = 2 , and 
magnification bias and evolution bias were neglected, 
s(z) = f^voi^) = 0 - The onset scale of the modifications 
in Eqs. (22) and (23) was chosen to be A = 10, which en¬ 
sures that they only become relevant at fc < 10“^ Mpc“^, 
i.e. beyond the matter-radiation equality peak in the 
matter power spectrum. This approximately marks the 

^ http://class-code.net 


maximum scale so far probed by galaxy surveys, as well 
as the transition away from the sub-horizon, quasi-static 
regime. One could of course choose a larger value of A 
to study effects on smaller scales, but we leave this to 
future work. 

Each of the panels in Fig. 4 shows the auto-spectra 
for four redshift bins, z = {0.1,1,2, 3}, all with tophat 
selection functions of full width Az = 0.2. The z = 
0.1 redshift bin consistently shows the smallest deviation 
from GR, which is straightforward to understand in light 
of Fig. 1 - at such small distances from the observer, 
even the largest angular scales correspond to subhorizon 
physical scales, where the MG corrections are suppressed 
in our parameterisation. The picture is more complicated 
at higher redshift however, as the relative importance of 
the various contributions to Eq. (28) varies non-trivially 
with scale and redshift. For example, the top right panel 
of Fig. 3 shows that the MG functions {fi, 7 ) grow with 
redshift for the GPL ansatz, and yet the z = 3 auto¬ 
spectrum deviates from GR less than the z = 2 curve. 
This is not the case for the fg\ = 0.1 curves, where the 
z = 3 curve is larger than the z = 2 curve for i < 10. 

This behaviour can be explained by inspecting the rela¬ 
tive size and differing modifications to the various terms 
in Eq. (28). The density term grows with time and is 
generally the dominant contribution to the total (see 
the left panel of Fig. 5), but this is not to say that its 
deviation from GR grows with time, or that it is the 
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Fig. 6.— Fractional deviation of the number count power spectrum for the A (left) and CPL-like (right) parameterisations for the z = 1 
and z = 3 redshift bins (indicated by light/heavy curves), with different values of the magnification bias (solid/dashed/dotted curves). The 
s = 0.4 (dashed) curves correspond to the intensity mapping case. 


dominant contribution to ACe- The right-hand panels of 
Fig. 5 show the deviation from GR of the A-ansatz power 
spectra, with the density-only and several non-density 
terms from Eq. (28) switched on and off. For ffiA = 0.1 
(left column) and z = 0.1 and 1, the total AC( curve 
(solid black) follows the behaviour of the density-only 
curve (solid blue). The non-density terms are generally 
subdominant to the density term (see the left panel), so 
while their fractional deviations can be large (e.g. see 
the leasing term; short dashed lines, right panels), their 
absolute contribution to ACe is small at these redshifts. 
This is not the case at z = 3, where the total devia¬ 
tion departs from the shape of the density deviation; the 
non-density terms are more important at this redshift, 
and the density-only deviation is close to zero, so does 
not contribute to the total deviation anyway. Note that 
cross-correlations between the density and non-density 
terms further complicate this analysis, as do differences 
between the two ansatzes (e.g. compare the very different 
behaviours of the density. Dm, in Fig. 3). Still, it is clear 
that the lensing and RSD terms are more important in 
driving the deviation from GR at high redshifts, making 
these interesting effects to focus on when searching for 
possible modifications to GR on ultra-large scales. Other 
terms also have significant fractional modifications, but 
are strongly subdominant for all but the lowest multi¬ 
poles; we study the modifications for the full set of terms 
in Eq. (28) in Appendix C. 

Another way of getting a handle on the relative im¬ 
portance of the various terms is to consider different 
values of the magnification bias, s(z), which multiplies 
some, but not all, of the terms. Fig. 6 shows the rela¬ 
tive deviations for both ansatzes at z = 1 and z = 3, 
for constant s = {0.0, 0.4,0.8}. Recall that s = 0.4 is 
the value to use when considering brightness tempera¬ 


ture fluctuations (e.g. for intensity mapping), for which 
lensing, ISW, and some other terms completely cancel. 
The magnification bias varies significantly between dif¬ 
ferent galaxy surveys (and also with redshift), but is less 
than 0.4 for most future surveys (Alonso et al. 2015). 

In Fig. 6, the z = 1 curves for = 0.1 and f^i = 0.1 
(lower panels) barely change with s, suggesting that the 
s-dependent terms contribute little to the total devia¬ 
tion. This is not the case at z = 3, where (e.g.) the 
integrated lensing contribution is expected to be more 
important, although the differences are still small. The 
contrast between the /^a = 0.1 and ffxi =0.1 curves 
(upper panels) is more striking, however. For the GPL 
ansatz, the dependence on s is very weak for both red¬ 
shift bins, but the A ansatz is strongly s-dependent. In 
fact, the z = 3 curve for /^a = O.I goes to zero when 
s = 0.4, suggesting that the deviation depends almost 
entirely on the lensing term. This is borne out by Fig. 5, 
which shows that the lensing term has a relatively large 
magnitude at z = 3, and sustains the largest fractional 
modification of any of the other (important) terms. 

4.3. Observability of deviations from GR 

The upshot of the above is that the distinctive de¬ 
pendence of ACi on different MG parameterisations and 
choices of MG ansatz, the physical parameters of the sur¬ 
vey (e.g. bias), and redshift, could in principle be used 
to differentiate between different MG theories. The devi¬ 
ations from GR are rather small for the MG parameters 
we have considered here, however, with a typical magni¬ 
tude of ^ 0.1% at £ < 10 for the A ansatz, and closer 
to 0.5% for the GPL ansatz. The deviations increase 
towards lower i for almost all redshift bins, typically be¬ 
coming important only for £ < 30, and they all but vanish 
beyond £ = 100. For reference, the cosmic variance error 
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Fig. 7. — Fractional deviation of the (unlensed) CMB tempera¬ 
ture power spectrum from its GR behaviour, for both parameter- 
isations. The results are comparable for the lensed temperature 
spectrum. Broadly speaking, the CPL ansatz shows larger de¬ 
viations as it affects higher redshifts than the A ansatz, thereby 
boosting the integrated Sachs-Wolfe effect. 

for a full-sky experiment in a single redshift bin is about 
30% at £ = 10, and 10% at £ = 100. The question, then, 
is whether it is even possible to observe deviations from 
GR at this level. 

It is first prudent to check whether the modifications 
considered above could have stronger effects on other ob¬ 
servables that might be easier to measure. The most ob¬ 
vious place to check is the CMB temperature (TT) power 
spectrum, which has been measured down to the cosmic 
variance limit for ^ < 1500 over almost the entire sky 
(Planck Collaboration XV 2014). Our parameterisation 
leaves early times unmodified by construction, but the 
observed TT spectrum still depends on late-time physics 
through secondary anisotropies like lensing and the ISW 
effect. Fig. 7 shows the deviation of the CMB TT spec¬ 
trum from GR for the same sets of MG parameters pre¬ 
sented in Fig. 4. The differences are again negligible for 
all but the lowest I modes, but are substantially larger 
than for the number count power spectra, weighing in 
at around 5-10% at £ < 10. Nevertheless, they are still 
well below the cosmic variance error, and so would not 
be detectable. 

While we only observe one CMB, it is possible 
to observe multiple redshift bins with multiple sur¬ 
veys/tracers, each of which can yield additional informa¬ 
tion about possible deviations from GR. It is therefore 
not clear whether the effects shown in Fig. 7 are too small 
to be observable due to the fundamental limit imposed 
by cosmic variance. To get a crude estimate of their 
detectability, we can write the modified number count 
power spectrum as 

= [1 +Agi{zi)]C^^(zi,Zj)6ij, (29) 

where we have explicitly neglected cross-correlations, and 
assume that the deviations can be parameterised by an 
amplitude parameter. A, and a known redshift and an¬ 
gular scaling, gi{z). We can roughly approximate the 
scaling as being constant with redshift and a step func¬ 
tion in £, such that the deviation is the same in each of 



Fig. 8. — Fractional deviation of the number count power spec¬ 
trum at 2 : = 3 for the A parameterisation (//IA = 0.1), as a function 
of the scale of the modifications, A. 

Vbins redshift bins for all £ < fmax- For a noise-free ex¬ 
periment, the single-parameter Fisher matrix for A can 
be approximated as 

FaA « ^/skyiVbins , (30) 

and the relative error on the amplitude is 

— ~ ^ -H \ 

A ^ Y AkyiVbinA^max + iy ■ ^ ^ 

For the CPL ansatz, the typical deviation is A ^ 0.3% 
out to £max 30. For a full-sky experiment (/sky = 1), 
one would therefore need > 200 redshift bins to obtain 
better than a 100% error on the amplitude parameter, 
which is clearly unrealistic (especially when one takes 
into account the fact that the size of the deviation will 
change for different redshift bin widths). 

This simple estimate does not give the full picture, as 
we have concentrated only on one set of survey charac¬ 
teristics, as well as neglecting noise and a number of nui¬ 
sance parameters. Certain types of survey might exhibit 
larger deviations, for example due to a different mag¬ 
nification or evolution bias. Since at least some of the 
terms in Eq. (28) depend on the particular properties of 
the tracer population (e.g. 6, s, /^o): it is also possible 
to use the ‘multi-tracer’ effect (Seljak 2009) to circum¬ 
vent the cosmic variance limit, as discussed in Lombriser 
et al. (2013). Multi-tracer analyses also appear to be 
necessary to detect other ultra-large scale effects, like the 
scale-dependent bias due to primordial non-Gaussianity 
(McDonald & Seljak 2009; Dalai et al. 2008; Ferramacho 
et al. 2014). Indeed, most of the relativistic terms in the 
number count spectra are undetectable in single-tracer 
analyses (Yoo et al. 2012; Alonso et al. 2015). 

We have also considered only a single set of MG pa¬ 
rameters, and certain viable theories might predict signif¬ 
icantly larger (or smaller) effects. We have argued that 
modifications around the Hubble scale (A ^ 1) with a 
similar redshift scaling to dark energy are in some sense 
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‘natural’ if the MG theory has been invoked to explain 
cosmic acceleration, so the choices we made in this re¬ 
gard should be somewhat representative. We did make 
relatively arbitrary choices of the amplitude parameters, 
/a .7 = 0-1) fixed A to a slightly higher value of 10, 
however, which may not be so typical. Fig. 8 shows the 
effect of changing A, but a more detailed investigation of 
the parameter space is necessary to draw more solid con¬ 
clusions (e.g. see Lombriser et al. (2013) for the effects 
in some example theories). Nevertheless, it seems likely 
that the deviations of the observables will be small, and 
thus difficult to detect, in most cases. 

5. CONCLUSIONS 

An impasse is emerging in the use of conventional 
galaxy surveys to test ideas about cosmic acceleration. 
As constraints from the linear, subhorizon regime are cy¬ 
cled back into MG models, theorists will likely begin to 
focus on theories whose phenomenology becomes man¬ 
ifest only outside the linear regime. Experiments are 
pushing towards high-accuracy measurements of smaller 
scales, but interpretation of the data there remains chal¬ 
lenging due to nonlinearity and poorly-understood bary- 
onic physics, even for ACDM. Tests of gravity on these 
scales rely heavily on computer simulations, and while 
rapid progress is being made (Winther & Ferreira 2014; 
Barreira et al. 2014; Falck et al. 2015), these remain 
resource-intensive and difficult to validate. 

Given this situation, it is logical to consider exploring 
the other end of the distance ladder. Although accu¬ 
rate measurements of horizon-scale modes are, experi¬ 
mentally, slightly further off than small-scale advances, 
they could ultimately prove more useful for shaping ideas 
about alternative gravity theories. If a theory’s effects 
are ruled out on horizon scales, then it becomes ques¬ 
tionable whether it can be of any use in driving cosmic 
acceleration, irrespective of small-scale phenomenology. 
An additional advantage is that, from the perspective of 
using effective field theory (EFT) to explore MG mod¬ 
els, large cosmological scales are safely in the low-energy, 
linear regime. This guarantees that one has a reliable 
theoretical understanding of the EFT, and concrete pre¬ 
dictions can be made. 

This is not to say that ultra large-scale cosmology will 
not face its own challenges. A particular obstacle is that 
any new horizon-scale physics emerging at late cosmo¬ 
logical times only imprints its effects on a small number 
of angular multipoles. As we have shown in this paper, 
deviations from the GR number count power spectra are 
;< 1% for ~ 10% corrections to the Poisson equation 
and slip relation, assuming that the MG effects occur on 


scales just above current observational reach. Deviations 
of this magnitude will be dominated by cosmic variance 
errors, although there is some hope that this could be 
ameliorated through the use of multi-tracer techniques. 
It is also possible that specific gravity theories may be 
found that behave very differently from the generalised 
parameterisation used here, resulting in more readily de¬ 
tectable deviations. 

The contribution of relativistic corrections to the angu¬ 
lar spectrum of number counts leads to interesting scale- 
and redshift-dependence of the MG deviations. Some of 
the terms (e.g. lensing) are also modified more strongly 
than others, which is potentially very useful if they can 
be separated out. It will be a struggle to isolate most 
of the relativistic corrections from the contribution of 
the matter density perturbation, however, even with the 
SKA (Maartens et al. 2013; Alonso et al. 2015). The 
best prospect is likely to be a multi-tracer analysis; the 
cross-correlation of an SKA 21cm intensity mapping sur¬ 
vey with photometric galaxy redshifts from LSST is suf¬ 
ficient for ~5 — 20cr detections of the amplitude of the 
lensing term and the combination of the other relativis¬ 
tic terms, for example (Alonso & Ferreira 2015). If we 
find in future that the relativistic effects are more clearly 
discernible than expected, then this could also be strong 
evidence that there is some large-scale modification to 
gravity boosting these terms. 

One could draw a parallel between our investigation 
of large-scale MG phenomenology and that of direct de¬ 
tection of (non-primordial) gravitational waves (GWs). 
Though an initial detection of GWs has yet to be made, 
the effects of modified gravity on compact object inspi¬ 
rals are investigated (e.g. Berti et al. 2013), and model- 
independent tools for GW analysis are constructed (Gor- 
nish et al. 2011) in preparation for a future era where 
GW detection is well-established. In a similar manner, 
this paper has presented an exploration of how simple 
changes to the linearised gravitational field equations 
propagate through to shifts in large-scale matter den¬ 
sity perturbations and related observables. As observa¬ 
tional constraints on scale-dependent growth continue to 
tighten, we are forced to ‘scan the horizon’ in this way 
if modified gravity is to remain a viable (but ultimately 
detectable) contender for explaining cosmic acceleration. 
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APPENDIX 

A. LARGE-SCALE DEVIATION OF 5Dm 

Here we provide a few more details on the linearised deviations from ACDM-I-GR on superhorizon scales that were 
discussed in §3.3. First we convert the independent variable to a; = In a, denoting derivatives with respect to x by a 
prime. Eq. (19) becomes 
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It is helpful to introduce the variable 5X = 5Y'yi/Y 'm ■ Making use of Eq. (18), we have 
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Integrating this in conformal time and substituting into Eq. (20), we then obtain 
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(A5) 


B. CONSERVATION OE 7^co 

In this appendix, we show that the comoving curvature perturbation is conserved on superhorizon scales in the 
presence of our modifications to the linearised field equations. The comoving curvature perturbation at late times, 
when pressureless matter is dominant, is 

7^co = -$ - jVm. (B1) 


Differentiating this and using the second of Eqs. (12) gives 

7t,o = -(4 + HT) + ^(H2-H). (B2) 

After some straightforward manipulation, one can eliminate Vm in favour of metric potentials. We differentiate our 
parameterised Poisson equation (1) and use the definition of Tm and the fluid conservation equations to obtain: 
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where C = d/dr] [ln(a^pM/i] ■ Rearranging, we obtain 
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Substituting this into Eq. (B2), we finally get 
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In the limit —>■ 0, the RHS of the above expression vanishes identically if we impose the reasonable restriction that 
fl does not tend to zero on ultra-large scales. The variable TZco is therefore conserved on superhorizon scales in our 
generalised scenario, in agreement with the findings of Wands et al. (2000) and Bertschinger (2006). 

Now that we know TZco —>■ 0 on large scales, Eq. (27) is straightforwardly obtained by substituting the —>• 0 limit 

of Eq. (B5) and its time derivative into (B2). 


C. TERM-BY-TERM DEVIATIONS FROM GR 

Figs. 9 and 10 show the fractional deviation of the number count power spectrum for individual terms in Eq. (28). 
Results are shown for two redshift bins, for each of the MG ansatzes considered above. The top panels show the 
absolute size of each of the terms in the GR case. We follow the naming convention for those terms given by Alonso 
et al. (2015); see Section 2.1.5. of that paper for the explicit expressions. 

Note that we only plot the autocorrelation for each term, e.g. the ‘Density’ curve shows the density-density part 
of the number count power spectrum only. We have neglected to show cross-terms here, even though (e.g.) the 
density-RSD term is known to be significant. 
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Fig. 9.— Top panels; Fractional contribution of various terms from Eq. (28) to the total angular power spectrum, C^", for the 2 = 1 (left 
column) and 2 = 3 (right column) redshift bins. Middle and lower panels: Relative deviation of each of the terms from its GR value, for the 
A-ansatz (middle) and CPL-like ansatz (lower), again for the 2 = 1 and 2 = 3 redshift bins. The solid black line shows the total deviation 
from GR, Note that we have taken the absolute value of the relative deviation (some of the deviations are negative). 


We briefly note three features of these figures. First, it is clear that the density term always dominates the total 
power spectrum at these redshifts, with the lensing and RSD terms being the only other significant contributions. All 
other terms are significantly smaller (except for at the very lowest £ ^ 5), and decay rapidly with £. 

The deviation of the density term from GR does not necessarily dominate the total deviation, however, as in some 
cases is very small. For example, for /^a = 0.1 at z = 3 (Fig. 9, top right deviation panel), the density term 

never deviates from GR by more than 5 x 10“^%. While the absolute magnitude of the lensing term is around an 
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Fig. 10.— The same as Fig. 9, but for the remaining set of terms from Eq. (28). The various spikes are due to ACe crossing zero. 

order of magnitude smaller than the density term, its fractional deviation is 3 orders of magnitude larger. It therefore 
dominates the total fractional deviation, as can be seen from the similar shapes of the ‘Leasing’ and ‘All’ curves. 

Finally, figure Fig. 10 shows that many of the subdominant ‘potential’ terms, while very small in absolute terms, are 
modified by the largest fractional amount. For f^A = 0.1 at z = 3, for example, the ‘Potential 3’ term (the smallest 
of all the terms, oc is modified by almost 100% at £ « 10. While this dramatic fractional modification is 

unobservable in this context, it may be worth looking for other observables that are more sensitive to the absolute 
value of this term. 
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D. MODIFICATIONS TO CLASS 


The angular power spectrum calculations in this paper were performed using a modified version of CLASS 2.4.2, a 
publicly available Boltzmann code (Bias et al. 2011). The code is written in such a way that only a few modifications 
to the source terms of the perturbation evolution equations are needed - everything else is sufficiently general to 
accommodate our MG corrections at late time. We introduced our changes into the Newtonian gauge part of the code 
only (CLASS can use both Newtonian and Synchronous gauge), and allowed our modifications to switch on only at 
late times, z < 6, falling back to the default CLASS code at higher redshifts. 

The metric potential $ is evolved using an equation that depends on the GR field equations by default. We replaced 
this with an expression derived by differentiating the modihed Poisson equation with respect to conformal time, and 
then substituting Eq. (12) to get 


$ = 



2)F 




-1 





Q 'lj2 


(Dl) 


which is valid at late times for the standard fluid stress-energy sources. The other metric potential, T, is derived 
from the solution for <I> using the slip relation, which we also replaced with the modified version. We validated our 
modifications by confirming that the code gives the same results as the unmodified CLASS code in the GR limit, 
/I = 7 = 1. We also wrote an independent Python implementation of the MG evolution equations, and conhrmed that 
the results agree with the evolution calculated by the modified CLASS code. 

The precision settings of the CLASS integrator must be increased considerably beyond their defaults to obtain 
sufficient accuracy for the power spectrum calculations (i.e. to obtain an accuracy much better that the typi¬ 
cal ^ 10“^ differences between GR and MG spectra found above). This can increase running time considerably, 
but is necessary to avoid the results being dominated by numerical noise. We increased the Bessel function sam¬ 
pling to selection_sampling_bessel=3.0, reduced the Fourier-space step size in the transfer function calculation 
to k_step_trans_scalars=0.3 to improve the resolution of the integrator, and set the tolerance of the perturbation 
solver to tol_perturb_integration=le-7. We found that these choices give sufficiently stable results that do not 
change appreciably with higher precision settings. We also set k_scalarjiiax_tauO_over_l_max = 4 and switched off 
the Limber approximation at all (. values, which improves the precision at higher 

Our modifications to CLASS are publicly available at https://gitlab.com/philbull/mgclass. 
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